/*
 * Integration.cpp
 *
 *  Created on: Aug 11, 2010
 *      Author: bernardo
 */

#include "Integration.hpp"

void Integration::readParameters(std::istream& inFile) {
	fileValues.clear();
	fileValues.addParam("v",v);
	fileValues.addParam("D",D);
	fileValues.addParam("dPhi",dPhi);
	fileValues.addParam("dT",dT);
	fileValues.addParam("finalT",finalT);
	fileValues.addParam("dTSave",dTSave);
	fileValues.loadParamVars(inFile);
	nStages=v.size();
	if (nStages!=D.size()) {
		std::cerr << "v.size()!=D.size() in Integration::prepare\n";
		exit(1);
	}
	fileValues.addParam("nStages", nStages);
}

void Integration::writeParameters(std::ostream& outFile) {
	fileValues.writeParamVars(outFile);
}

void Integration::prepareDynamics() {
	std::vector<StageParam> stageParams(nStages);
		for (size_t i=0;i!=nStages;i++) {
			stageParams[i].setV(v[i]);
			stageParams[i].setD(D[i]);
		}

	dynamics.setStages(nStages,1/dPhi);
	dynamics.setStageParams(stageParams);
	dynamics.setdT(dT);
	double prob[] = {1,0};
	dynamics.density.deltaFunction(0.1+dynamics.density.getDX()/2,prob);
}

void Integration::evolution() {
	prepareDynamics();
	if (densFile!=NULL) {
		writeParameters(*densFile);
		dynamics.density.saveProbDens(*densFile,dynamics.getTime());
	}

	if (statFile!=NULL) {
		writeParameters(*statFile);
		dynamics.density.saveStat(*statFile);
	}

	int nSaved = 0;
	double finished=0;
	double outflow = 0;
	while (dynamics.getTime()<=finalT) {
		if (statFile!=NULL) {
			dynamics.density.calcStats();
			dynamics.density.saveStat(*statFile, dynamics.getTime());
		}
		if (densFile!=NULL && nSaved*dTSave<=dynamics.getTime()) {
			std::cout << dynamics.getTime() << "\t" << finished << "\t" <<
					1-finished << "\t" << dynamics.density.prob() << std::endl;
			nSaved++;
			*densFile << "\n";
			dynamics.density.saveProbDens(*densFile,dynamics.getTime());
		}

		//outflow = dynamics.oneStepUpWind(0);
		outflow = dynamics.oneStepCrNic(0);
		finished += dT*outflow;
	}
}

void Integration::completionTime(ProbDens<1>& compTime) {
	prepareDynamics();
	compTime.clear(0.01);
	double completed=0;
	while (1-completed>1e-10) {
		double outflow = dynamics.oneStepUpWind(0);
		completed += dT*outflow;
		if (completed>1e-10)
			compTime.add(dynamics.getTime()-dT/2,outflow*dT);
	}

/*	if (compTimeFile!=NULL) {
		fileValues.writeParamVars(*compTimeFile);
		completionTime.calcStats();
		*compTimeFile << "#prob: " << completionTime.prob() << "\n";
		*compTimeFile << "#meanT: " << completionTime.meanX() << "\n";
		*compTimeFile << "#varT: " << completionTime.varX() << "\n";
		completionTime.saveProbDens(*compTimeFile,0);
	}*/


}





